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The dynamical equations for periodic systems under constant external stress (E-print, arXiv:cond- 
mat/0209372 ) have been extended to many-body potential case. Rather than introducing La- 
grangian dynamics, Newtonian dynamics was employed directly. The kinetic energy term of stress 
and the imagined force only due to particles' passing through geometric planes or boundaries were 
also discussed. 
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The molecular dynamics (MD) simulation [1-3] is an important tool and widely used in research in Condensed 
Matter Physics, Biophysics, Chemical Physics, Space Physics, and even Astrophysics. The primary idea is to simply 
apply Newton's Second Law onto every particle. For simplicity or feasibility, the periodic boundary condition is often 
imposed, which results in that the system is modelized as a crystal in each simulation time step. Then only the 
MD particles (in the MD cell) and the periods are the degrees of freedom. The principle physics problem becomes 
dynamical equations for the degrees of freedom with the possibility of external stress/forces acting on surfaces of the 
system. Obviously Newton's Second Law is absolutely applicable directly onto the MD particles. What about the 
periods? 

By extending Andersen's idea [4] and introducing an unproved Lagrangian, Parrinello and Rahman proposed their 
theory regarding the dynamical equations for the periods in 1980 [5,6]. Although their theory has been extensively 
employed in simulations, we recently found that the MD particles violate Newton's Second Law in their theory [7]. 
The reason is that the kinetic energy for the MD particles in their Lagrangian is incomplete. Why such a Lagrangian 
Q ■ is adopted? 

The whole story is closely associated with another key problem: whether the internal stress contains the so called 
kinetic energy term. As essentially stress is the net force per unit area, it definitely contains the regular force term. 
_^ ' The truth is that the kinetic energy term is widely believed and calculated as the other part of stress. Irving and 
, Kirkwood proved the kinetic energy term for stress with pure statistics in 1950 [8]. Haile proved it by introducing an 
— i ' imaged force appearing only if a particle passing through a geometric plane in 1992 [3]. Zhou recently proved that 
the kinetic energy term was not physically grounded and should be deleted from the stress expression [9] . We made 
some discussion on physics of the imaged force in the TRANSPORT OF MOMENTUM section in our previous work 
[10] and found that for a system defined with a fixed space there should be a momentum flux term to cancel the 
unneeded momentum change per unit time only due to particles' passing through the system boundary in order to 
have Newton's Second Law applied correctly [10,11]. If we put this momentum flux term onto the momentum change 
side of the dynamical equation, nothing new; if we move it onto the force side, it can be regarded as an effective 
force and further results in the kinetic energy term of stress in principle. When we regard fluids as composed of 
independent cells (no periodicity) and separate their whole motion into macroscopic and internal ones, the kinetic 
energy term of stress can be explained as momentum flux due to particles' internal motion across boundaries or 
collisions of internal motion between particles from a statistical point of view [11], resulting the stress similar to that 
proposed by Irving and Kirkwood [8]. The periodic system simulated with molecular dynamics as a crystal can be 
regarded as a special fluid, in which the periods represent the macroscopic motion and the motion of the MD particles 
is the internal motion. Then the kinetic energy term of stress from the MD particles can be suggested. Parrinello 
and Rahman obtained a stress with a similar kinetic energy term by choosing the above incomplete Lagrangian [5,6], 
but as a result, the MD particles violate Newton's Second Law [7]. If the complete Lagrangian were used there, the 
similar kinetic energy term of the stress would have not been generated in physics, as in our previous works [12,13]. 

Instead of introducing a Lagrangian, we applied Newton's Second Law onto the MD particles directly avoiding the 
above problem and applied Newton's Second Law onto halves of the system for dynamical equations of the periods 
(a, b, and c), and got the kinetic energy term of the internal stress by considering the collisions between particles and 
the boundary recently [10]. However this work was done for pair-potential only. Now let us extend it to many-body 
potential cases. The many-body force and potential were already discussed previously [13,14] with some misprint in 
Reference 13, but starting from a Lagrangian for constant-pressure molecular dynamics. Here we will start from the 
forces for the constant-stress molecular dynamics. For simplicity, we use Fig. 1 in version 6 of Reference 10 as Fig. 1 
here. 

For up to M-body interactions considered, the kernel task for this extension is to find the net force acting on the 
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half-line-cell bar B^ by all particles in the L part (left half) of the crystal in Fig. 1 

M 

F h = £ F h m) > (1) 

m=2 

where F^" 1 ' is the contribution of m-body interaction, and any period h = a, b, or c. For any particle Ik with 
k = 1, 2, • • • , or M in the whole crystal, its position vector is vi k — T/ fc + r^, where is the position vector for its 
image particle in the MD cell, and 

Ti k = / fe , a a + 4, b b + 7 fe , c c (2) 

is the lattice translation vector of the cell in which it resides. For m-body forces, the interaction between the right 
and left half of the crystal means all participating particles must be distributed in both halves, or all participating 
particles must not be in any one of the halves. Then is the summation of all these possibilities devided by the 
number of equivalent bars to Bh in the right half crystal iVh = |S| / |<7h|, where S is the cross section area vector 
in Fig. 1, and Oh = dfl/dh with cell volume ft. Now let us write out the force for all the possibilities with total t 
particles (m > t > 1 ) distributed in the right half of the crystal and the rest particles in the left half 

Cl,h,/2,h,-,/t,h>0) (/t+l,h,i"t+2,h,---,-fm,h<0) t 

F S=^77^tV E E £ f r(^ M ,---,r,J, (3) 

y '' i ./..•••./•; {i t+1 ,i t+2 ,-,i m } ,u=i 

where we use to denote the summation over indexes in {• • •} with conditions in (■■■), and 

f j™^ (r/j , r j 2 , r j 3 , • • ■ , r j to ) is the m-body force acting on particle I M . We also require that any two of the indexes 
in {• • •} must not refer to the same one physical particle, while any particle and its image particle in any other cell 
are different physical particles. Remembering 

/l,h>0 /i.h>0 

E — E EE E - w 

II i"i,h Ii,h' h,h" il£Cell M D 

with h', h"= a, b, or c, but h' ^ h, h" ^ h, and h" ^ h', iVh = J2i , J2i ,, 1> an ^ the crystal translatability, we 
have 

. +OC (/2,h,/3,h,"'Jt,h>0) (i"t + l,h,J"t+2,h,---,J"m,h<0) t 

FS-^TT^vE E E E f rK' r ^^,--,r,J (5) 

' V '' /i,h=0 {ii,I 2 ,h,-,It} {h+i,h+2,-,I m } M=l 

1 +OC (I2,hj3,h,---Jt,h>0) (h+l,h,h + 2,h,---Jm,h<0) 

E E t f iT ) ( r h^i2^hr--,ri m ). (6) 
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With the crystal translatability and moving the system so that the cell containing particle I\ is translated to the MD 
cell, it becomes 



— OO (^2,h,^3,h,---,^t,h>0 (^t+l,hi-ft+2,hr"!-fm,h<i) 

F *"^ = tUm^W E E E r\rur l2 ,r l3 ,...,r Im ), (7) 

' v >' l=0 {i,I 2 ,I 3 ,-,I t } {I t +l,I t +2,-,I m } 



and Eq. (I) becomes 

M m— 1 



F * = E E F S (8) 

m=2 i=l 

M m— 1 — OO (^2,h|J3,h>"')^t,h>0 (It + l,\i,It + 2,\i-----Jm,h<l) 

= E E a _ lV(m _ fV E E E f- (m) (r 2 ,r /2 ,r /3 ,-..,r /m ). (9) 

TO=2 t= l V ^-V /• ( =0 {i,I 2 ,I 3 ,-,I t } {I t+1 ,I t+2 ,-,I m } 

Meanwhile the total up to M-body cell potential energy is 
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where the m-body cell potential energy 



E ™ = 7J E ^ m) (r 4 ,r /2 ,r 73 ,...,r 7m ), (11) 



m! r 

{«,7 2 ,7 3 ,--.,7 m } 



with m-body potential (^ m )(r 7l , rj 2 , • • • , rj m ). This form of E 1 ^™/ ensures that only a fraction ^ of </?( m ' belong to 
the cell when k (m > k > 0) particles are in the cell. Now let us take a derivative 



-k 1 ®- . £ w£w*.**,-,*.). (i2) 



d „f m i 1 

{ij2,7 3 ,---,7 m } 

with 

f j fc m) (r 7l , r Ja , • ■ ■ , r Jm ) = - V r , fc ^ (m) K , r 7a , • ■ ■ , r Jm )■ (13) 
The right side of Eq. (12) can be split into two terms based on the sign of 7 m> h 
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where 

7 m , h >o 
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F K=^rW E Wi: ) (r i ,r, 1 r„-,r / J (15) 

^ +00 An,h>Z 

= m(m-2)' E E fff^.^^/s.-.^J. (I 6 ) 

V ; ;=0 {i,7 2 ,7 3 , -,I m } 



7 m , h <0 

F ^=^TW E Wr^r^r,^...,^) (17) 

V h {i,I 2 ,I 3 ,-,I m } 



_^ — CO Im , h < 1 



TWE E flr^r,,^,^,-..,^). (18) 



m(m- 2)! 

v ; ;=0 {i,7 2 ,7 3 ,---,7 m } 

Eq. (16) can be expanded regarding t' , the number of particles distributed in the right half of the crystal, of particles 
indexed from I 2 to I m -1 



m-2 +oo 7 m . h >i (-f2,h,73,h,"-,7 t / + lh >0) (l t i +2 ,h A'+3,h,'"Jm-l,h<0) f (mw 

F S= E £ £ % *L7-w ■ < 19 » 

t'=0 (=0 7 m {i,7 2 ,7 3 ,-,7 t , + 1 #7 m } {7 t , +2 ,7 t , +3 ,...,7 m _ 1 } 



Then let us use the translatability, move the system so that the cell containing particle I m is translated to the MD 
cell, and rename particle I m as particle i and particle i as particle I m 



m-2 -oo ( 7 2,h,7 3 ,h,---,7 t / + lih >;) (7 t / +2 ,hi J t'+3 ,hi---. 7 m-i,h,7 m , h <() (m) / 

„(m) 1 v-^ tj ('Ti,r l2 ,r l3 ,- ■■ ,r Im ) 

*m-- m 2^2^ 2^ 2^ t'\(m-t'-2)\ [M) 
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Please note the m-body force and potential are independent on the order of the m particles' position variables in the 
expression of f j™^ (r j x , r j 2 , • ■ • , r j m ) and </?( m ' (r/ x , r j 2 , • • • , r j m ) respectively. 

For each given I (negative) value, Eq. (18) can be expanded with t' particles of those indexed from I 2 to l m -\ 
distributed in the right half of the crystal defined by Ik,h > I 



m-2 -OO ( 7 2,h,J"3,h,--->J"t' + l,h>0 ( J "t'+2,h> J "t'+3,h>--- Jm-l,h,-f m ,h<l) „ ( m ) , . 

p(m) _ ^1 \ ^ J m ^ I r/2 ' r/3 ' ' ' ' ' Vlm I 

h - ~ 7n ^ ^ ^ ^ V\(m-t'-2)\ ' 

'- 0,=0 {i,l2,h,-J t , + 1 } {/ t / +2 ,/ t / + s,-,/m-l,/m} 

Making use of crystal translatability again, and 



(22) 



^^(r/^r/^r^.-^r/J = 
for conservation laws of energy and momentum, we have (I < 0) 

(l2,h,h,h,---,It' + l,h> 1 ) ( J l'+2,hi'('+3,h)" , i'm,h<l) 

E E -(m-t' -l)f^ ) (r i ,r l2 ,r l3 ,---,r Im ) 

{ 7 t'+2. / t'+3.'-' Jm} 

E E (t' + l)tt l \^r l2 ,r l3 ,---,r Im ), 

{i,l2,h,---,I t '+i} { 7 t'+2^t'+3'---. / ™} 



then 



As a result, we have 
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F h= E E - E ( F K + *t>) = -^c«. 



m— 2 i— 1 m— 2 

With the main interaction tensor defined for up to M-body cases 
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The m-body contribution can also be written as 

— (m) 1 
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The general interaction tensor should be 
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r= e ^ (m) . (32) 

m=2 



where the m-body contribution can be extended as 

->-(m) 1 ™ 

£ m = V^\ E S f ir ) ( rj i' r ^' r/ 3'---' r ^) (n M -ro) (33) 

' {/iJ 2 ,/ 3 ,---J m }A'=l 

= Vm\ E ^ f i^\ r h,ri a ,T Ia ,---,T Im )T Ill (34) 
' {h j 2 , i 3 ,- ,i m }n=i 

for an arbitrary given position vector ro, with volume V of the whole system and Eq. (23). With translatability, it 
can then be written as 



,(m) 



{ij2,-"J m }M=2 

{i,/2,-,/m} 



+ i E fi m) (r„r /2 ,...,r 7m )r 4 + 

{»,/2,-,/m} 
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+ E E ^ " T,„, r, 2 - T,„ , • • ■ , r M • ■ • , r Jra - T,, )r„ 

_>._>.(m) _>_>.(m) 

£ m + £ p . (35) 



where with n, the total number of particles in the MD cell, 



1=1 

Fl m) = -^"T'^^^^'Vj = -V r ,l£g. (37) 

^ {/2,/3,-,-fm} 

By applying Newton's Second Law onto the right half of the crystal and applying statistics over all possible collisions 
between particles and boundaries and over all possible cut bulks of complete cells from the same crystal shown in Fig. 
2 of version 6 of Reference 10, the same form of dynamical equations are generated for many-body interactions as in 
Section V. of version 6 of Reference 10. For an example, Eq. (56) of version 6 of Reference 10 is still 

WH = (^T + "iQ • ~cT, (38) 
where W is the equivalent mass for periods, T is the external stress, H = {a, b, c}, a = {er a , (7b, <r c }, and 

7T =£m + £ p +^E Wl l^' 2 ' ( 39 ) 
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with particle mass and e p = J2m=2 £ p For ^ ree systems, as T = and no collision, Eq. (38) becomes 

vph = ^ e„T +i J2 F ^ ■ ~^ ( 4 °) 

where 

M 



F i= ^ F,[ m) = -V ri £ ceH . (41) 



m=2 

As stated above, the MD particles obey Newton's Second Law 

m i r i =F i (i = 1,2,- ■■, n) . (42) 
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